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Abstract. We study the propagation of an unusual type of periodic travelling 
waves in chains of identical beads interacting via Hertz's contact forces. Each 
bead periodically undergoes a compression phase followed by a free flight, due to 
special properties of Hertzian interactions (fully nonlinear under compression and 

Ph ■ vanishing in the absence of contact). We prove the existence of such waves close to 

binary oscillations, and numerically continue these solutions when their wavelength 
is increased. In the long wave limit, we observe their convergence towards shock 
profiles consisting of small compression regions close to solitary waves, alternating 
with large domains of free flight where bead velocities are small. We give formal 
arguments to justify this asymptotic behaviour, using a matching technique and 
previous results concerning solitary wave solutions. The numerical finding of such 
waves implies the existence of compactons, i.e. compactly supported compression 
waves propagating at a constant velocity, depending on the amplitude and width 

f—^ | of the wave. The beads are stationary and separated by equal gaps outside the 

wave, and each bead reached by the wave is shifted by a finite distance during a 
finite time interval. Below a critical wavenumber, we observe fast instabilities of 
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rvq ■ the periodic travelling waves leading to a disordered regime. 
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1. Introduction and main results 



Understanding wave propagation in granular media is a fundamental issue in 
many contexts, e.g. to design shock absorbers [HI HE], derive multiple impact laws 
[25| [3J [351 [331 El], detect buried objects by acoustic sensing [41)] . or understand 
possible dynamical mechanisms of earthquake triggering [32J. One of the important 
factors that influence wave propagation is the nature of elastic interactions between 
grains. According to Hertz's theory, the repulsive force / between two identical 
initially tangent spherical beads compressed with a small relative displacement 5 is 
f(S) = k5 a at leading order in 5, where k depends on the ball radius and material 
properties and a = 3/2 (see figure [1]). This result remains valid for much more 
general geometries (smooth non-conforming surfaces), and a can be even larger in 
the presence of irregular contacts [311 [22]. The Hertz contact force has several 
properties that make the analytical study of wave propagation more difficult than 
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in classical systems of interacting particles : the dependency of f(S) on 5 ~ is 
fully nonlinear for a > 1, /"(0) is not defined for a < 2, and no force is present 
when beads are not in contact. This makes the use of perturbative methods rather 
delicate, since the latter often rely on nonlinear modulation of linear waves which 
are not present in such systems, and usually require higher regularity of nonlinear 
terms. 




Figure 1. Schematic representation of two identical and initially 
tangent spherical beads that are compressed and slightly flatten, the 
distance between their centers decreasing by 5 ~ 0. 

The simplest model in which these difficulties show up consists of a line of identical 
spherical beads, in contact with their neighbours at a single point when the chain is 
at rest. For an infinite chain of beads, the dynamical equations read in dimensionless 
form 

— y = V'(x n+ i - x n ) - V'(x n - x n _i), n e Z, (1) 

where x n (t) e K. is the displacement of the nth bead from a reference position and 
the interaction potential V corresponds to Hertz contact forces. It takes the form 

V(x) = -^—\x\ 1+a H(-x), (2) 

1 + a 

where H denotes the Heaviside function vanishing on M~ and equal to unity on IR + 
and a > 1 is a fixed constant. System ([2]) is Hamiltonian with total energy 

^ = £^) 2 + ^ + l-Xn). (3) 

For a = 3/2, Nesterenko analyzed the propagation of compression pulses in this 
system using a formal continuum limit and found approximate soliton solutions with 
compact support J37J EH] (see also jUHUJIH] for more recent results and references). 
As shown by MacKay [36] (see also [30]), exact solitary wave solutions of ([I]) exist 
since an existence theorem of Friesecke and Wattis [21] can be applied to the chain 
of beads with Hertz contact forces (see also [IB] for an alternate proof). Moreover 
these solitary waves have in fact a doubly-exponential decay [HI US, SB] that was 
approximated by a compact support in Nesterenko's analysis. 

Much more analytical results are available on nonlinear waves in the granular 
chain when an external load /o is applied at both ends and all beads undergo a 
small compression 5q when the system is at rest. Around this new equilibrium 
state, the dynamical equations correspond to a Fermi-Pasta-Ulam (FPU) lattice 
[23J [39], that sustains in particular stable solitary waves with exponential decay 
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(see e.g. (2TJ [191 123 [20], [26] ) and periodic travelling waves also known as "nonlinear 
phonons" [T7J [271 1131 ESS 121] • The existence and qualitative properties of periodic 
travelling waves in granular chains are important elements in the understanding of 
energy propagation and dispersive shocks in these systems, as shown in reference 
[TJ] in the similar context of FPU lattices. However, the persistence of periodic 
travelling waves when f — > is not obvious because the sound velocity (i.e. the 
maximal velocity of linear waves) vanishes as f , and the uncompressed chain of 
beads is commonly denoted as a "sonic vacuum" [38J. 

In this paper we show that periodic travelling waves with unusual properties 
exist in granular chains under Hertz contacts without precompression. These waves 
consist of packets of compressed beads alternating with packets of uninteracting 
ones, so that any given bead periodically switches between a free flight regime and 
a contact regime where it interacts nonlinearly with two or one neighbours. Waves 
of this type have been numerically computed in reference J57J , for periodic granular 
chains consisting of three or four beads. Here we provide an existence theorem 
valid at wavenumbers close to tt and proceed by numerical continuation when the 
wavenumber goes from it to the long wave limit q — > 0. The limit q = tc corresponds 
to binary oscillations, i.e. travelling waves with spatial period 2 where nearest 
neighbours oscillate out of phase. In the long wave limit, the periodic waves display 
small compression regions close to solitary waves, alternating with large domains of 
free flight where bead velocities are small. 

Our existence result for periodic travelling waves close to binary oscillations is 
described in the following theorem. In what follows, Cp er (0, 2n) denotes the classical 
Banach space of 27r-periodic and C k functions u : R — > R, endowed with the usual 
supremum norm taking into account all derivatives of u up to order k. 

Theorem 1. System (QP with interaction potential (0|) admits two-parameter fami- 
lies of periodic travelling wave solutions 

x^(t;a,fj) = aw jU [(7T + /j,)n±a~ T ~t], (4) 

parametrized by a > and \x G V, where V denotes an open interval containing 
0. The function u^ is odd, 2-n -periodic and satisfies the advance- delay differential 
equation 

<(£) = VM£ + vr + /i) - «„(£)) - HMO - «/.(£ - * - aO). ( 5 ) 

The map fi i— > u^ belongs to C^V, Cp er (0, 2tt)) and satisfies the symmetry property 
u n(0 = ~ u -n(£, + 7r )- The function Uq is determined by the initial value problem 

< + W'(u ) = 0, (6) 






«o(0) = 0, u' (0)= Po , Po = (l + a)^ w 7 " , (7) 

1+cJ 

where W(x) = Y^\ x \ 1+a is a symmetrized Hertz potential and T(x) = f Q °°e' t t x ~' 1 dt 
denotes Euler's Gamma function. Moreover, for fi < 0, u^ is a linear function of £ 
on an interval [— £iOt/),£i(ju)] with £i(/-0 = M/2 + °(/ i )- H takes the form 

MO = P^for all £ e [-Zi(jt),Si(ji,)], (8) 
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with p,j, = p + 0(\fj,\). When £ = (n + /i) n ± a 2 £ g (— £i(aO + 2/c7r, £i(/i) + 2&7r) 
('A; e Zj, t/ie nt/i particle performs a free flight characterized by 



a + l 



x n+1 (t;a,/x) > x n (t;a,/i), x n (£; a,//) > x n _i(t; a,//), x n (t; a,fi) = ±a 2 p M . 

To prove theorem [Q we locally solve equation (jSj) using the implicit function 
theorem, where we have to pay a particular attention to regularity issues due to the 
limited smoothness of the interaction potential V. This approach is reminiscent of 
the previous work [28], where we proved the existence of periodic travelling waves in 
the Newton's cradle, a mechanical system in which the granular chain is modified by 
attaching each bead to a local linear pendulum. However, in that case the periodic 
travelling waves were obtained by nonlinear modulation of linear oscillations in the 
local potentials, a limit which does not exist in the present situation. 

A complementary approach to the analysis of exact periodic travelling waves 
consists of obtaining approximate solutions described by continuum models. Ap- 
proximate periodic travelling waves have been previously studied through a formal 
continuum limit in the granular chain (see J3H], sections 1.3-5, [IE] , section 4.6 and 
references therein), but in a case where all beads were interacting under compression. 
Closer to our case, Whitham's modulation equation may be used to approximate 
the solutions of theorem [1] and study their stability properties, following the lines of 
reference [13] where the example of a granular chain with linear contact interactions 
was analyzed. 

A limitation of theorem [1] stems from the fact that it provides an existence result 
for wavenumbers q = n + fj, close to 7r. To analyze the existence of periodic travelling 
waves on a full range of wavenumbers, we numerically continue (for a = 3/2) the 
solutions of theorem [1] by decreasing the wavenumber q down to values close to 0. 
We compute the nonlinear dispersion relation associated to periodic travelling waves, 
and analyze their asymptotic form in the long wave limit q ~ 0. In this limit and 
when the wave velocity is normalized to unity, one finds large ensembles of beads 
performing a free flight at small velocity, separated by smaller regions consisting 
of 4-5 balls under compression with relative displacements close to a solitary wave. 
Using a matching technique, we perform a formal asymptotic study of the limit 
q — >■ which provides a limiting wave profile in good agreement with our numerical 
computations. 

These results are completed by a stability analysis. We compute the Floquet 
spectrum modulo shifts of the periodic travelling waves, and find fast linear in- 
stabilities below a critical wavenumber q c ps 0.9. This threshold corresponds to 
the case when the average number of interacting beads in the compression regions 
becomes larger or equal to three, or equivalently to the maximal number of adja- 
cent interacting beads becoming > 4. In this regime, the initial numerical errors 
made on the travelling wave profiles are amplified during dynamical simulations, 
leading to the disappearance of the travelling waves after some transcient time. A 
disordered regime is established shortly after the instability, but interestingly some 
partial order can be still observed in the form of intermittent large-scale organized 
structures. In contrast to the above situation, long-time dynamical simulations of 
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periodic travelling waves with wavenumbers q > q c yield wave profiles that remain 
practically unchanged when propagating along the lattice. The existence of very 
slow instabilities in this regime is a more delicate question which will be addressed 
in a future work. 

Besides periodic waves, our numerical results also demonstrate the existence of 
compactons (solitary waves with compact supports) in chains of beads separated by 
equal gaps outside the compression wave. For q < q^ w 1.8, the periodic solutions of 
(J|)J) obtained numerically behave linearly on intervals of length larger than the delay 
q = 7T + jjl involved in (j5]), and these solutions can be linearly extended in order to 
get new solutions of (jSJ). Using this property and the Galilean invariance of ([I]), 
we obtain compacton solutions for which the beads are stationary and separated 
by equal gaps outside the wave, and each bead reached by the wave is shifted by a 
finite distance during a finite time interval. 

Our situation is quite different to what occurs in the absence of gaps between 
beads, because in that case compactons exist in continuum models of granular chains 
[37J [38l H], k^ a transition from compactons to noncompact (super-exponentially 
localized) solutions occurs when passing from the continuum model to the discrete 
lattice [3]. In our case the existence of compactons is linked with the occurence of 
free flight, made possible by the unilateral character of Hertzian interactions which 
vanish when beads are not in contact. 

The outline of the paper is as follows. In section [2] we prove theorem [1] and study 
some qualitative properties of the travelling wave profiles. The numerical results are 
presented in section [3j which contains in addition a formal asymptotic study of the 
long wave limit and the discussion of the existence of compactons. Section @] gives 
a summary of the main results and mentions implications for other works. Lastly, 
some useful properties of solitary wave solutions are recalled in an appendix. 

2. Local continuation of periodic travelling waves 
Periodic travelling wave solutions of ([1]) take the form 

&n(*) =«(£)» £ = qn-ut, (9) 

where u is 27r-periodic, q G [0, 2n) denotes the wavenumber, u G M \ {0} the wave 
frequency and (el the spatial coordinate in a frame moving with the wave. Thanks 
to a scale invariance of ((Tj), the full set of periodic travelling waves can be deduced 
from its restriction to u = 1. Indeed, due to the special form of the Hertz potential 
(J2J), any solution x n of flTJ generates two families of solutions x[ if) = ax n (± a^~ t) 
parametrized by a > 0. Consequently, any solution of the form (^ with u = 1 
corresponds to two families of periodic travelling waves propagating in opposite 
directions 

a>)(£) = au(qnTa^t), (10) 



a-l 



with frequency to = ± a a ' 

Fixing u = 1, equation ([1]) yields the advance-delay differential equation 

«"(0 = V'(u(Z + q)- «(£)) - V'(u(Z) - «(£ - q)), £ G R, (11) 
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with periodic boundary conditions 

w(£ + 2tt) =«(£)■ (12) 

Equation (TTTj) possesses the symmetry u(£) — > — w(— £) originating from the reflec- 
tional and time- reversal symmetries of (JTJ). In the sequel we restrict our attention 
to solutions of (TTT1) invariant under this symmetry (i.e. odd in £). This assump- 
tion simplifies our continuation procedure, because it eliminates a degeneracy of 
periodic travelling waves linked with the invariance of (fTTj) under translations and 
phase-shift. 

2.1. Binary oscillations. The following lemma exhibits a particular solution of 
(1 11 1) - 0121) obtained for q = it and corresponding to binary oscillations in system ([]]). 

Lemma 1. For q = it, there exists an odd solution uq of < f77]) -< f7^j determined by 
the initial value problem (E))-([?p, and satisfying 

uo(£ + 7r) = -u (£). (13) 

In addition, there exists a one-parameter family of solutions of (TJP taking the form 

x n {t) = a(-l) n+1 w (a^£), aG(0,+oo). (14) 

Proof. We look for solutions of ( ITT]) satisfying ffTB"]) . so that the advance-delay dif- 
ferential equation reduces to the ordinary differential equation (J6J) with W(x) = 
^V(—2\x\) = Y^\%\ 1+a - Now we must find a 27r-periodic solution uq of <Q satisfy- 
ing (TT3|) . Equation fl6]) is integrable since / = \{u') 2 + W{u) is constant along any 
solution u, and its phase-space is filled by periodic orbits. Any solution u of ([6]) 
with initial condition w(0) = 0, u'(0) = p > is odd in £ (due to the evenness of W) 
and its period T(p) is given by 

/"^max 

T{p)=4 {p 2 -2W{v))- 1/2 dv 



Jo 
where u max = |[p 2 (l + «)]*+«. Our aim is to find p > such that T(p ) = 2n. 

1 — q 1 — q 

Using the change of variable v = u max t 1 +°' , one finds T(p) = 2(1 + aji+a f3p 1 + a ) 
where 

1 L 



/? = f t^- 1 {l-t)- 1/2 dt = B( 
Jo 



l + a'2' 



and B{z 1 w) = yu+w) denotes Euler's Beta function (see [2], formula 6.2.1 p. 258). 
Since T(l/2) = y/ir, one obtains finally 

T(p) = cap&, c a = 2y^(l + a)^ r _T~° lv (15) 

Consequently, for p varying in (0, +oo) the period T(p) is strictly decreasing from 
+oo to 0, and there exists a unique po > such that T(p ) = lix. For the corre- 
sponding solution uq of (jSJ), ^o(O) = po is given by ((7j). Lastly, there remains to 
check property ( ITB1 in order to ensure that uo satisfies (TTTj) . Using the fact that 
m (vt) = (since uq is odd and 27r-periodic), we have u' q {tt) = —p by conservation 
of /. Consequently, m (£ + tt) an d «o( — are solutions of the same Cauchy problem 
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at £ = 0, which implies that f lT3|) is satisfied. Solutions ( }T4|) are obtained using 

(Ho]). □ 

Remark 1. The solution uq of lemmaU\is non-unique because —uq is also a solution, 
but this second choice corresponds to shifting solution JT^j j by one lattice site (or 
equivalently performing a half-period time-shift). In addition, the solutions of fifty 
with period 2n/{2k + l) (k G N) also satisfy fnj\) , and thus they are solutions of < f77)j - 
( TJJ| ) for q = 7T. However, these additional solutions do not yield any new travelling 
wave because they can be recovered from u$ by tuning a in (fl^| ). Moreover, the 
solutions of IHty- fnfy with period ir/k do not satisfy fiTBj). hence they do not satisfy 

(Cnp-mp /or ?=7r. 

2.2. Periodic travelling waves close to binary oscillations. In order to locally 
continue the solution u of (fTT|) -( lT2"|) for q m n we need to define a suitable functional 
setting. In the sequel we note X k = Cp er (0, 27r) and consider the closed subspace of 

X k = {ue C p fc cr (0, 2tt), w(-0 = -«(0 }• 
We denote by r q G £(Xfc) the shift operator (r q u)(£) = «(£ + g) (note that r^ maps 
Xfc into itself). Problem ffTl]) - ffT2|) restricted to odd solutions can be rewritten 

u" + N(u, q) = in X , u G X 2 , (16) 

where the map N(., q) defined by 

N(u, q) = V'((I - r_» - V'((r q - I)u) 

maps Xq into itself. The smoothness of N restricted to X2 x 1 is proved in the 
following lemma. 

Lemma 2. The map N belongs to C l (X 2 X M.,Xq) and 

D u N(u , tx) = i W"(u ) (I - r ff ). (17) 

Proof. Since V' G C 1 (IR), the map it h-> V'(u) belongs to C 1 (X ) (this classical result 
follows from the uniform continuity of V" on compact intervals). Consequently, to 
prove that N is C 1 it suffices to show that the map G : X 2 x M. — > X defined by 
G(u,q) = r q u is C 1 . Since the continuous bilinear mapping IT : X 2 X £(X 2 ,X ) — > 
X , (u,A) h-> Am is C°°, it suffices to check that the map g i-> r g belongs to 
C^R, £(X 2 , X )). Taylor's formula yields 

.. T q+h U - T q U ,., \h\ 

|| 7 T qU \\l°° S ~7T~ || It ||l°° 

for all u G X 2 and g,/ieR with /i ^ 0. Consequently, lim ~ (r g+ /j — r g ) = r g -fL in 

ft— )-0 * 

£(X 2 ,X ), i.e. -p- = T q j-; in £(X 2 ,X ). Moreover, by the mean value inequality 



\\i~ q +hu' - T q u'\\ L °o < \h\ ||n"|| L oo, 
hence ||t 0+ j, 4? — r a -4 II < \h\ and -?■ is continuous in the operator norm. 

This completes the proof that N is C . Formula ( ITTjl follows from elementary 
computations, using the chain rule, the fact that r n = r—n on X^ and t^uq = —uq, 
and the equality W"{x) = 2(V r "(2z) + V"{-2x)). D 
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In order to apply the implicit function theorem to equation ([TBI) in a neighbour- 
hood of the solution (u, q) = (u , ir), we now prove the invertibility of the operator 
L = £ + D v N(uo,ir). 

Lemma 3. The linear operator L G £(X 2 ,X ) is invertible. 

Proof. For a given / G X , we look for y G X 2 satisfying 

^1/ = /- (18) 

We use the splitting X k = X k © X^ , where 

X± = {ueX k , u{£ + tt) = ±u(0 } 

and denote by P 1 * 1 = |(J ± r w ) the corresponding projectors on X k . Since t*o G X^ 
and W is even, it follows that W'(w ) G X + and D u N(u , vr) = W"(«o) P~ G 
£(X ,Xq~). Setting f± = P ± f, y± = P ± y, problem f|T8|) can be rewritten 

y'i = U, (19) 

L-y- = /-, (20) 

where L~ y_ = y" +W"(u ) y_. We first observe that the operator jL G C(X^, Xq) 
is invertible, with inverse A G C{Xq,X^) given by 

O4/)(0 = - F {s-ir)f{s)ds+ f {t-8)f(8)ds (21) 

Wo Jo 

(injectivity comes from the fact that functions in X k are odd, and it is lengthy but 
straightforward to check that (I2TJ) satisfies the required periodicity and symmetry 
properties when / G Xq). Consequently, equation ( 1T9|) admits a unique solution 
y + G Xj". We now consider the case of equation (12"U1) . The operator T G £(X^", X^) 
defined by Ty_ = W"(uo) 2/- is compact (since X^ is compactly embedded in Xq~), 
hence L~ G C(X^ , Xq) is a compact perturbation of an invertible operator. As a 
consequence, L~ is a Fredholm operator with index and its invertibility will follow 
from the fact that Ker L~ = {0}. The proof of the injectivity of L~ is classical (see 
e.g. [16] p. 693) but we include it for completeness. In the sequel we denote by 
v p the solution of (jSJ) with initial condition v p (0) = 0, v'(0) = p, so that u = v Po , 
v p is odd and periodic with period T(p) given by ( 1T5|) . In addition we rescale v p 
in a 27r-periodic function w(r,p) = v p (-^-r). The solutions of the homogeneous 
equation y" + W"(uo)y = are spanned by two linearly independent solutions 

u' n and z = -^ , where u' n is even, 27r-periodic and z is odd. Since we have 
dp |p=po' o ' F 

v p (£) = w(^p- £,p) and T(p ) = 2tt, it follows that 

*(0 = —rfa) e<(o + ^(e,Po) (22) 

27T ap 

with T'(p ) = n (1 Z2) • Since the nondegeneracy condition T'(p ) =£ is satisfied, 2 
is the sum of an unbounded function of £ and a 27r-periodic one. As a conclusion, 
Ker L~ = {0} since 2 is non-periodic and u' is even, hence L is invertible. □ 

Lemma [3] allows one to solve equation (I16p locally in the neighbourhood of (u, q) = 
(uq, 7r) using the implicit function theorem. This yields the following result. 
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Lemma 4. There exists an open neighbourhood Q of uq in X 2 , an open interval I 
containing q = tt and a map $ G C l (X, X 2 ) (with $(71) = u ) such that for all gel, 
equation l[Tb}) admits a unique solution u in Q given by u = $(g). 

The function u^ describing the family of periodic travelling waves of theorem [1] is 
defined as 

u M = <l>(7r + /i). (23) 

Now there remains to analyze in more details the qualitative properties of u^. 

2.3. Qualitative properties. We start by examining a symmetry property of ffTTl) . 
One can check that u is an odd solution of (fTTj) - ([T2"l) with q = TT—fi if and only if — r w u 
is an odd solution for q = tt + yU. Consequently, by lemma H] (which ensures the local 
uniqueness of the solution) one has ^(tt+h) = —T n ^(TT—/j), i.e. m m (0 = —U-^+tt). 

The next lemma describes the monotonicity properties of u M in different intervals. 

Lemma 5. For all \x small enough, there exists 6(fi) G (0, 7r) such that u' (£) > 
for £ G [0, 6(h)), <(0 < for £ G (0(/x), tt]. Moreover lim 9(ji) = tt/2. 

Proof. From the construction of w in lemma [TJ the results holds true for \x = with 
9(0) = it/2, and we have u'q = —W'(uq) < on (0, 71"). Since w M is C 2 -close to w for 
/i ~ 0, for all e G (0, 7i"/2) one has u"^ < on [| — e, | + e] for /i small enough, m^ > 
on [0, f — e] and u^ < on [| + e, 7r]. Then the sign properties of u^ stated in lemma 
|5] follow from the monotonicity of u'^ on [| — e, | + e] and the intermediate value 
theorem, and the asymptotic behaviour of 8(/S) is obtained by letting e — > 0. D 

Since u^ is odd and 27r-periodic, it satisfies w M (0) = m^(vt) = and u^(tt + h) = 
—u^(tt — h). By lemma it follows that u^ is positive on [0, tt], maximal at C, = 9(fi) 
and negative on [tt, 2tt]. Figure [2] illustrates the profile of m m for fi w —0.63. 

Using lemma [5j we now prove the existence of an interval where the function u^ 
is linear. 

Lemma 6. For all \x < small enough, there exist £i(/x) G (0, #(//)) snc/i t/iat 

^(e + vr + /i)-^(0 >0 /oraM£e(-7T-^-£i,6), (24) 

«,*(£ + *• + /*)-«/*(£) <0 for all £e {&,<*- H-&), (25) 

ei(A*) = -f + o(A*)- (26) 

In addition, one has 

MZ)=P^forall£e[-Z 1 ,£ 1 ], (27) 

wt/i p^ = p + 0(|/i|) and p defined by (0). 

Proof. We first show that equation 

« M (6 + vr + /i) = «„(&) (28) 

admits a solution £i(/i) G (0, 0(h)), a property clearly illustrated by figure |2J From 
lemma |5] it follows that for all £i G [0,0], the equation 

« M (0=<V(£i) (29) 
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Figure 2. Graph of u^ for q = tt + // rj 2.51, numerically computed 
from equation (ITT1) using the method described in section [3j The 
profile of u^ is linear on the interval [— £i,£i] (see lemma EJ), and the 
values of w M (£) at £i and £ 2 = £i + Q coincide. 



admits a unique solution £ = £ 2 in [9, tt], the latter depending smoothly on £i G [0, 0) 
by the implicit function theorem. Moreover, d(£i) = £ 2 — £1 is a strictly decreasing 
positive function of £1, with d(0) = n and d{6) = (see figure EJ). Consequently, for 
/i G (— 7r,0), the equation 

d(£i) = tt + // (30) 

admits a unique solution £i(/x) in (0, #), which depends smoothly on ji by the implicit 
function theorem. Moreover, expansion (126]) is obtained by solving ( 130]) for (£ 1; fj,) rj 
(0, 0) with the implicit function theorem, and using the identity 



rf'(0) 



u '^) 



-l = -2 + 0(|/i|), 



(31) 



where we note p M = w^(0). The first equality in ( 131]) is obtained by considering the 
case £1 ~ of equation ( 129|) . and the second one holds true because u' (tt) = — u' (0). 
By definition of d, £i(//) defines a solution of (1281) and £1 (/i) + 7r + /i G (6 1 , tt). 

In the sequel we prove inequalities ( 124)) and (125]) . The proof is based on the 
known monotonicity properties of u^, and can be followed more easily with the help 
of figure |2j 

Let us first prove (EIJ by treating the cases £ G (— g — £1, — £1] and £ G (— £i,£i) 
separately (we recall that q = tt + \x G (0, 7r)). 

Choosing £ G (— g — £1, — £1] yields — £1 < £ + q < 2n — q — £1, therefore we 
have u M (£ + g) > inf[_£ Lj27r - (? -£ 1 ] u M = w M (-£i). Since w M (-£i) > w M (£) for all £ G 
(-9 - £1, -£1], we have m^(£ + g) > u M (£). 
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In the second case £ G (— £1, £1), we have £1 < £ + g < £1 + g for /x small enough, 
hence w M (£ + q) > inf [ ?li§1+g] m m = w M (£i). It follows that u M (£ + q) - m m (£) > 
«„(&) - «„(£) > for all£ e (-6,6)- 

Consequently, we have proved inequality (T2"lj) . Now let us prove (1251) , treating the 
cases £ G (£i, £i + q] and £ G (£1 + q, 2n — q — £]_) separately. 

For £ G (£i, £i + g] we have £i + g<£ + g<27r (since £i + q < n and q < 7r), hence 
M (£ + q) < su P[^ 1+q ,2n] u » = <v(£i + q)- Consequently we have u M (£ + q) - m m (£) < 
uJii +q)~ m m (0 < for £ G (£i, £i + q], hence w M (£ + q) - w M (£) < 0. 

In the case £ G (£i + g, 27r — q — £i), we have 2n — £i — q < £ + g < 2n — £i for 
/i small enough, hence u M (£ + q) < sup [27r _ 6 _ (/)27r _ ?i ] u^ = u^tt - £i - q). Since 
m m (2tt - £i - g) < ?v(£) for £ G (£i + g, 2vr - g - £x), we have m m (£ + g) < w M (£), 
which completes the proof of inequality ( ]2"5|) . 

Now we can deduce (J2"T|) from the advance-delay equation (1TT1) and inequality 
( |24|) . The latter implies m m (£) —u^ — q) > for all £ G (— £i,£i + g). Consequently, 
for all £ G (— £i, £i) we have both « M (£ + g) — w M (£) > and m^(£) — w M (£ — g) > 0. 
Equation (II ip yields in that case w^(£) = 0, from which (12TJ) follows easily. D 

As a conclusion, from lemmas HJ |6] and expression ({TO]) , we have obtained the 
two-parameter families of periodic travelling wave solutions x^(t; a, /i) of theorem [I] 
close to the binary oscillations described in lemma HJ 

Now let us discuss in more detail some implications of lemma |6j In what follows 
we shall write x^(t;a, /i) = x n (t) for notational simplicity. Since x n (t) = aM M (£) 

a — 1 

with £ = (tc + ft) n±a^~t, it follows that as soon as £ G (— £i(/i) + 2A;7r,£ 1 ( / u) -\-2kii) 
(k G Z) we have 

x n+1 (t) > x n (t), x n (t) > x n -i(t), x n (t) = ±a~ p M . (32) 

In other words, each bead periodically switches between a free flight regime corre- 
sponding to ( 1321) and a contact regime where it interacts with two or one neigh- 
bours. From inequalities (12^1) and (1251) . interaction with two neighbours takes place 
for £ G [£i + 7T + fi + 2kTT, -k — \i — £i + 2/ctt]. This completes the proof of theorem HJ 

Remark 2. 5y Galilean invariance of (QJ), we deduce another family of solutions 

x n (t) =x n (t) +vt, v = Ta SL 2~p^ (33) 

In that case, each bead periodically switches between a pinning regime where it re- 
mains stationary, and a contact regime where it interacts with two or one neighbours. 

1 — a 

Each bead is shifted by =)=27rap M after one cycle of period T = 2na~^~ , therefore 
bead displacements are unbounded in time. Such wave profiles are reminiscent of 
stick-slip oscillations occuring in the Burridge-Knopoff model [13] , albeit this model 
corresponds to completely different mathematical and physical settings involving dif- 
ferential inclusions and nonlinear friction. 
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3. Numerical results 

The analysis of section [2] has proved the existence of periodic travelling wave 
solutions of ([1]) with wavenumbers q ~ it, close to the binary oscillations ( flU) cor- 
responding to q = tt. In this section we numerically continue this solution branch 
by decreasing q up to values close to and analyze the qualitative properties of the 
waves. Our computations are performed for a = 3/2. In section I3"U| we deduce the 
existence of compactons from the numerical results obtained for periodic waves. 

3.1. Numerical continuation. We discretize problem ( riT]) - ([T2"j) using a second- 
order finite difference scheme with step h = 7r.l0~ 3 , and consider the discrete set 
of wavenumbers q = m j^ obtained with all integers m in the interval [1, M], where 
we fix M = 50. With this choice, the corresponding delays q occuring in (flTT) are 
multiples of the step h, hence the advance-delay terms of ( ITTT) can be computed by 
a discrete shift of the numerical solution. In addition it is important to keep h <^i q 
due to a boundary layer effect that occurs at small wavenumbers (see below). The 
resulting nonlinear equation is solved iteratively using a Broyden method [12j and 
path-following, starting the continuation at q = tt and u = Uq. Note that higher- 
order finite difference schemes are useless for a = 3/2, because V G C 1 ' 1 ^ 2 (R) has 
only a limited smoothness, which guarantees that u is C 3 but not C 4 . 

Figure displays the solution u computed numerically for several values of q. As 
shown by figure IU the supremum norm of u diverges as q — > with 

2 

IMIoo ~ k a q^-« (34) 

and k 3 / 2 ~ 1.36. 

The solution is found exactly linear (more precisely, both terms at the right side 
of (fni) vanish) for f G \-£{q) + 2kn,i(q) + 2fcvr] (k G Z), with lining) = tt and 

q— >-0 

\im£(q) = 0. For q ps 0, the linear part of u behaves like w(£) ~ — q 1 ^ £ for 

q— >tt w 

£ G [—£(q),£(q)] (see figure [5]). More generally, the monotonicity properties of u 
proved in lemmas |5] and [6] for q = it + /i w it are still valid for the numerical solution 
with q G (0, it] (with the correspondence £(q) — £i(g — it) between the notations of 
lemma |6] and the present ones) . 

Since x n (t) = u(qn — t), the linear behaviour of u corresponds to a free flight 
of some beads, separated by equal gaps of size u'(0)q. Free flight occurs when 
qn-t G {-£{q) + 2kir,£(q) + 2kn) (k G Z), i.e. we have 

qn-t G (-£(q) + 2kir,£(q) + 2kiT) <^> x n+1 {t) > £„(£), x n {t) > x„_i(t), x n {t) = 0. 

(35) 
For £ G (0, 2tt), the two linear parts of «(£) are connected by an inner solution that 
becomes steeper as q — > 0, and extends over an interval of length 2(tt — £{q))- This 
part of w(£) corresponds to the nonlinear interaction of some packets of beads. 

It is interesting to study the size of the packets of interacting beads and beads 
in free flight as a function of the wavenumber q. Following (l35|) . and recalling that 
q = m -g, we define the fcth packet of beads in free flight as the set of beads with 
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index n in the interval 

hk(t) = (q-'i-iiq) +t) + — , q~ l {l(q) +t) + — ), (36) 

m m 

and the feth packet of interacting beads as the set of beads with index n in the 

interval 

hk+i(t) = [q-\i(q) +*) + —> <T M(?) + *) + ^ "]• ( 37 ) 

Note that Ik(t) moves with the wave, so that any given bead periodically switches 
between the free flight and contact regimes. 

Let us define Ofe(t) = q~H + ^ fe + q~ l d(q) and 

iV(g) = 2g- 1 (7r-%)), (38) 

so that l2fc+i(£) — [ a fc(^), Ofc(t)H-iV(g)]. We denote by jV(o, L) the number of integers 
in the interval [a, a + L], equal to [-^J or |_£J + 1 depending on the values of a and 
L (\_L\ denotes the integer part of L). The function M satisfies 

Af(a + l,L)=Af(a,L), / M(a,L)da = L. 

Jo 

The number of beads with index n in I 2 fc+i(£) is then equal to Af(ak(t), N(q)), and 

switches between the two values [N(q)\ and [N(q)\ + 1 during the motion. This 

number being g-periodic in t and m-periodic in k, the average number of adjacent 

interacting beads reads 

— i^f N(a k (t),N(q))dt= [ Af(a,N(q))da = N(q) (39) 

q m t^ 1 Jo Jo 

(we have used the changes of variables a = a,k{t) in above computation). Figure 
[6] depicts the evolution of N(q) when q is varying in (0,7r). One observes that 
N(n) = 2 (as expected for binary oscillations, since beads are paired for almost all 
times), and the value of iV at our minimal q is iV(7r/50) ~ 4.8, which is close to the 
spatial extent of Nesterenko's solitary wave [38] . 

In addition, using similar arguments as above (just considering l2k(t) instead of 
^2k+i(t))i t ne average number of adjacent beads in free flight is F(q) = — — , i.e. it is 
equal to the length of the interval of linear increase of u divided by the wavenumber 
q. The number of beads with index n in the open interval lakit) switches between 
the two values |~-F(g)] — 1 and \F(q)~\ during the motion, where \F] denotes the 
smallest integer not less than F (note that \F~\ — 1 and [-^J coincide when F is not 
an integer). Figure [7] provides the graph of F(q). One can notice that F(q) > 1 when 
q < q k ?z 1.8. In that case, two packets of interacting beads are always separated by 
some beads in free flight. 

Now let us describe the nonlinear dispersion relation satisfied by the periodic 
travelling waves. Let us denote by u(£; q) the family of 27r-periodic wave profiles 
parametrized by q G (0, 7r], computed by numerical continuation from u(.;tt) = u . 
We recall that these solutions correspond to travelling wave solutions of ([T]) given 
by 

x n (t) — au(qn — a~~ s ~t + <fi; q), (40) 
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Travelling wave profile 




Travelling wave profile 




Figure 3. Solution of (TTTj) - (fT2l) computed numerically for q 

(top), q = 37r/10 (middle) and q = 7r/50 (bottom). 



TV 



where a > is a parameter. Instead of a we now choose the wave amplitude 

A = ||{x n }||L<x,( ZxR ) = a \\u(.; 



as a new parameter. The frequency of solution ( HOI) is consequently given by the 
dispersion relation 



u(q,A) = A 3 *-\\u(.;q)\\ > 



(41) 
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Figure 4. Graph of the supremum norm of u (in logarithmic scale) 
when the wavenumber q varies. The line corresponds to the approxi- 



mation II u I 



k a g 1 -" for a = 3/2 and k 3 / 2 ~ 1.36 




Figure 5. The curve gives for the different wavenumbers q the 
constant value of u' in the free flight interval £ G [—£(q),£(q)} (plot in 
logarithmic scale). The line corresponds to the approximation u' ~ 






for a = 3/2 and k 3 / 2 ~ 1.36. 



From the above numerical results, we deduce 

Lo(q,A)r^k a A 2 q as q — > 0. 



(42) 



Consequently, the dispersion relation ( I4ip behaves linearly in q at fixed A^O when 
q — > 0, and one has |r(0, 0) = 0, which corresponds to a vanishing "sound velocity" 
in the limit of small amplitudes. This property is consistent with the fact that, 
in the granular chain with a precompression /o, the sound velocity vanishes when 



/o - > 0. The graph of u(.,A) is shown in figure M for A 

u(n,A ) = l). 



A 



\Uo\ 



(so that 
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Figure 6. Average number N(q) of consecutive interacting beads 
when the wavenumber q varies. 
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Figure 7. Graph of F(q) = — — , the average number of adjacent 
beads in free flight. 



In what follows we formally analyze why the scaling (|3"4"|) of \\u(.; q)^ occurs 
when q — > and heuristically compute the constant k a present in the dispersion 
relation 



3.2. Long wave limit. Let us examine the limit q — > more closely. We first 
normalize the solution of ( 1TT|) -( TT2|) computed in section |3~T1 by setting u = g 1 ^ u, 
which yields 



q 2 u"(0 = V'(u(£ + q)~ «(£)) - V'(u(0 - u(£ - q)), 



or equivalently 



u»(t) = q*' 1 - (V'[ ^ + q) ~ m ]-V'[ H0 - ^ - q) ]) , 
q q ' q 



(43) 



(44) 
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Figure 8. Plot of the dispersion relation (jUJ) for a = 3/2 and 
A = 1 1 Mo I |„o) an d its linear approximation at q = given by 



where we have used the fact that V'(x) 
right side of ([HI) by its continuum limit 



\x\ a H(—x). For q m 0, we replace the 



u"(0 = q a - 1 ^v\u'). 



(45) 



Setting q = in ( l4"5"j) and using the fact that u(0) = 0, we get the outer approximate 
solution 

u o (0 = z^ Ce[0,7T), (46) 



TV 



where A will be subsequently determined by a matching condition. When q is small, 
equation (1431) with g = does not describe the travelling wave profile in a boundary 
layer around £ = 7r where w' becomes large, and a rescaling of w is necessary to 
capture the structure of the inner solution. Setting s = q~ x (£ — 71") and £t(£) = y(s), 
equation (j4"51) becomes 

(47) 
lim y so i(s) = fc a , lim y so i(s) = -k a , (48) 

s— >— oo s— S-+00 



y"{s) = V'(y(s + 1) - y(s)) - V'(y(s) - y(s - 1)), s e 
There exists an odd solution y so \ of (14T1) satisfying 

fc a , lim y so i(s) = -k a , 

where the convergence towards fc a is super-exponential and 



k 



3/2 



1.3567 



(49) 



(see appendix [A]). This solution corresponds to a solitary wave solution of ([TJ 
[36| [15| 148] with velocity equal to unity. It yields the inner solution 



«t(0 = l/sol ( 



£-7U 



(50) 



which agrees very well with the numerical solution in the vicinity of £ = 7r (see figure 
|9]). Now let us match the outer and inner solutions at some point £ = tc — q A(q) 
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when q — > 0, where A(q) > 0, lim ? _ > o (qA(q)) = and lim^o A(q) = +oo. The 
matching condition reads 

limA(l-^(g)) = limy sol (-A(g)), 

q— >0 7T q— >0 

which gives consequently A = k a . This approximation is in reasonable agreement 
with the numerical results since the value of ks/2 deduced from figure |U differs from 
the value ( 14"9"1) by 8.10 -3 (q should be further decreased to get a better agreement). 
From the previous analysis, we deduce the following approximate solution of ( T4"3"j) 
obtained by summing the inner and outer solutions and substracting their common 
value at the matching point 




*w»(0«< I q P-« ( 51 ) 

Q 

The numerical and approximate solutions are compared in figure [10J which shows a 
good agreement between both profiles. The mathematical justification of this formal 
asymptotic analysis is left as an interesting open problem. 

With approximation (I5TJ) at hand, let us now examine the wave structure in more 

2 

details. Setting a = q^ 1 , = n in (j4"0]) and using (I5TI) yields approximate travelling 
wave solutions x^ v {t) = w ap p(g(n — t) + it) with velocity equal to unity and 0(1) 
amplitude when q — y 0. These solutions take the form 

< PP (t) = —q(n-t) + y sol (n - t), for \n - t\ < - (52) 

7r q 

and are —-periodic with respect to the moving frame coordinate s = n — t. These 
waves consist of a succession of compression solitary waves, separated by large re- 
gions of free flight (of size 0(q~ 1 )) where particles move at a constant 0(q) velocity 
and equal gaps of size 0(q) are present between adjacent beads. 

3.3. Compacton solutions. Several approximations of compression solitary waves 
with compact supports have been derived in the literature [37J [381 H] > m order to ap- 
proximate exact solitary wave solutions of ([!]) decaying super-exponentially [T51H8] . 
However, the existence of exact solitary wave solutions of ([1]) with compact support 
remained an open question. Such solutions are supported by different classes of 
Hamiltonian PDE with fully-nonlinear dispersion, and have been known as com- 
pactons after the work of Rosenau and Hyman [41] . Nonlinear PDE supporting 
compactons can be introduced to analyze lattices with fully nonlinear interactions 
near a continuum limit. In this context, a common scenario is the transition from a 
compacton to a noncompact (super-exponentially localized) solution when passing 
from the continuum model to the discrete lattice [1] (see also [121 [29] for similar 
results in the context of discrete breather solutions). 

In contrast to this situation, we show in this section that the numerical results of 
sections 13.11 imply the existence of compactons for the full lattice ([[]) . This result is 
due to the occurence of free flight analyzed prevously, made possible by the unilateral 
character of Hertzian interactions, which vanish when beads are not in contact. 
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Figure 9. Solution of ( l43j) with period 2tt computed numerically 
for q = 7r/50, with a zoom in the vicinity of £ = it (black curve). This 
solution is compared with the inner solution ( 1501) corresponding to 
the solitary wave studied in the appendix (blue curve). Both profiles 
agree very well in a domain of width 5 ~ 0.2 around £ = n, 8 being 
approximately four times larger than 1/q. 




Figure 10. Solution of (14311 with period 2ir computed numerically 
for q = 7r/50 (black curve) and its approximation (1511) (blue curve). 



Let us consider the solutions w(£; q) of (1TT]) - (fl2"]) obtained numerically. These 
solutions behave linearly on the intervals I = [—£(q),£(q)} and I\ = [2tt — £(q), lit + 
£(q)]- Denoting by p(q) = u'(0;q) their slope in these intervals, we have u((,;q) = 
p(q) £ on I x and u(f ; q) = p(q) (f - 2tt) on I 2 . 

In what follows we assume F(q) = 2£(q)/q > 1, which corresponds to fixing 
<? < Qk ~ 1-8 (see figure [7]). This condition can be interpreted as follows. In the 
case of a strict inequality F(q) > 1, in a chain of beads where the wave w(£; q) 
propagates, we have seen that two packets of interacting beads defined by (1371) are 
always separated by some beads in free flight, hence w(£; q) can be interpreted as 
a periodic train of independent pulses of finite width (i.e. compactons). Moreover, 
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under the assumption F(q) > 1, the length of the intervals I\ and I2 is larger than 
the delay q involved in (ITTj) . 

If F(q) > 1, the solution u(£; q) can be linearly extended in order to get a new 
solution of ( ITT]) , defined by 



p(q)Z 


for 


£ < -%), 


«(£;?) 


for 


eGH(g),2vr + %)] 


p(?) (f - 2vr) 


for 


£>2vr + %). 



^(£) = < «(£; g) for £ e H(g), 2vr + %)], (53) 

(p(q)(Z-2n) for £>2vr + %). 

The profile of U q is plotted in figure [TT] for q w 0.88. Let us check that C7 g defines 
a solution of fTTTj) when F(q) > 1 (for notational simplicity we shall omit the q- 
dependency). For £ G I~ = (-00, £] we have U"(g) = 0, U(£ + q) - U(£) > and 
U(£) — U(£ — q) > 0, hence C/ is a solution of (II ip on I~ . The same properties hold 
true for £ G J + = [27r — £, +00). Moreover, [/ is a solution of (ITT]) for £ G J c = 
[9 — £, 27r + £ — g] because U = u and r± g C7 = r± q u on J c . Since the condition F > 1 
is equivalent to having K = J - U J c U J + , U defines a solution of (11 ip on R. 

The solutions of ( ITT]) defined by (T55|) correspond to travelling wave solutions of 
(|T]) given by 

X n (t) = aU q (qn- a^t + (f)), (54) 

where a > 0, q < q% and G K are parameters. Moreover, by Galilean invariance of 
([TJ, we deduce another family of solutions 

X n {t)=X n {t)+vt (55) 

where we fix 

v = a 2 p{q)- (56) 

These solutions correspond to single compactons. Indeed, they consist of travelling 
waves, in the sense that bead velocities X n (t) and relative displacements X n {t) — 
X n _i(i) are functions of n — ct with c = q~ x a~^~ . Moreover, each bead remains 
stationary except in a finite time interval where it experiences a compression. For 
a bead with index n, compression occurs when qn — a^~t + <p G (£(q), 2tc — £(q)), 
i.e. when U q does not behave linearly. Fixing <fi = —£(q), the support of the moving 
compacton at time t corresponds to n G (ct,ct + N(q)), where N(q) is defined by 

Beads outside the support are separated by equal gaps of size A = ap(q) q and 
each bead position is shifted by 2iiap(q) after the passage of the compacton. The 
gap A depends on the parameters a and q, or equivalently on the compacton width 
N(q) and velocity c. Fixing c = 1 and letting q — > 0, the gap becomes 0(q) and the 
compacton profile converges towards the classical solitary wave in the compression 
region (see section 13~2"]) . 

Similarly to (153]) . one can obtain a family of 2-compacton solutions of (TTT]) by 
gluing the solutions U q and f/ g6 ;(£) = p(q) 9 + U q (£ — 2n — 9). Let us define 

r(2) f U q (0 for Z<2n + 9 + £(q), 



U q , e (0 H f/W^) for £> 27 r-^ (5?) 
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where 9 = qd, d > is a parameter and q G (0,5k], so that q < 2£(q). The graph 
of U q J is illustrated by figure [TT] (note that £/q and C/^ j coincide on [2n — £{q), 2tc + 

+ ^(5)])- Let us check that t/j^ defines a solution of (TTTj) (for simplicity we shall 
omit dependency in q and 9 in notations). Firstly, U^ is a solution of ( 111]) for 
£ G J~ = (-oo,2tt + 9 + £ - 5] because £/ (2) = C/ and r ±q U i2) = r± q U on J~. 
Similarly, U {2) is a solution of (JTTJ) on J+ = [2ir - £ + q, +00) because f/ (2) = £/ (1) 
and T± q UW = r± q U^ on J+. If 5 < £ + 5/2 we have J~ U J+ = R and thus 
f^ 2 ) defines a solution of (ITT|) on R. Moreover, if £ + 5/2 < q < 2£, then we have 
(C/( 2 ))"(£) = 0, f/( 2 )(e + 5) - f/ (2) (0 > and U {2 \0 - U^ 2 \^ - q) > for all 

1 E J c = (2tt + 9 + £ - 5, 2vr - £ + 5), hence f/ (2) is a solution of (ED on J c , which 
yields a solution of (III]) on R = J" U J c U J+. 



Each solution of (HI 



X n (t) = af/ (2 J(5n-aVt + 0)+ t ,t (5J 



with f defined by (156]) consists of two independent compactons separated by the 
distance ^(5) + d, with stationary beads outside the compactons. Similarly, gen- 
eralizing the above construction would allow to define an infinity of solutions of 
(PQ), corresponding to an arbitrary number of independent compactons with variable 
compacton spacing. 

Note that the remark made here on the existence of compactons follows from the 
numerical results obtained on problem (|TTl) - (fl2l) for wavenumbers q < q^. Obtaining 
an analytical proof of the existence of compactons remains an open problem. 

3.4. Error minimization. In this section we evaluate and improve the numerical 
precision of the computations performed in section 13.11 Our numerical approach 
is in the same spirit as former computational methods developed to analyze the 
propagation of discrete breathers in nonlinear lattices [5] (see also [IT]). 

We consider a chain of 100 particles with periodic boundary conditions 

x n+N (t)=x n (t), AT =100, (59) 

so that the lattice period is an integer multiple of the wavelength A = 27r/g = 100/m 
for our discrete set of values of q. We fix a = \\u(.\ 5)||~ in (HOI , where u(.; q) is the 
solution of ffTT]) - (fT2]) that was computed numerically for different values of 5. We use 
the resulting profile as an initial condition in system ([T]), which determines a periodic 
travelling wave solution X n = (x n ,x n ). Numerical integrations are performed using 
the standard ODE solver of the software package Scilab, where we have fixed the 
relative and absolute error tolerances to 10~ 12 and 10~ 14 respectively. For each 
solution u{.\ 5), we compute the relative residual error 

||{x,*.iCn-*»(o)}JL 



E(q) 



ll{*»(0)}J 



with T = qa 2 ( so that an exact travelling wave solution with velocity T x would 
cancel E(q). The inverse wave velocity T(q) is consequently given by 

T(q) = q\\u(-q)\\J 1 . (60) 
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Figure 11. Top panel : graph of a solution of ( ITTj) defined by (15311 . 
corresponding to a single compacton (we have set q ~ 0.88). Bottom 
panel : graph of a solution of ( ITT]) defined by (157)) . corresponding to 
a 2-compacton (we have fixed the same value of g and 9 = 2). 



Figure [12] displays the graph of E(q) (dotted line). With our numerical solution, 
this error remains less than h 2 ps 10 -5 for q > 67r/25 and grows when q is further 
decreased, reaching 3.10 -4 at q = 7r/50. We attribute these larger errors to the 
sharpness of the travelling wave in the compression region when q is sufficiently 
decreased. One way to improve the results would consist in discretizing ( TTT]) - (TT2"]) 
with a nonuniform mesh, thinner in the compression region. However, in what 
follows we use a different approach which decreases the residual error by several 
orders of magnitude on the whole range of wavenumbers q. 
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We use the Gauss-Newton method [12] in order to determine a refined initial 
condition {X n (0)} n=1 N minimizing ||{X ra+1 (T) — X n (0)} n || 2 (the Newton method 
is not directly applicable due to the noninvertibility of the Jacobian matrix at the 
exact solution). The improved initial condition obtained in this way for a given 
wavenumber q will be denoted by X n (Q) = X„ ' . To compute these solutions, we 
initialize the Gauss-Newton iteration using the numerical solutions u(.;q) of (ITTj) - 
(JT2J1 obtained previously, since they already provide initial conditions quite close 
to the optima. To reduce the computational cost, we modify the Gauss-Newton 
method by actualizing the Jacobian matrix only at some steps of the iteration, 
when a re-evaluation is required to decrease the residual error. The residual error 
E(q) obtained with this numerical method drops to 3.10" 10 or less, as shown in 
figure IT21 (full line). At the end of the Gauss- Newton iteration (Arth iteration), the 

||{x< fc )(o)-x< fc - 1 )(o)} n 



relative difference between the last two iterates r(q) 
below 4.10~ 7 (see figure fl2l dash-dot line). 



\\{xk h \o)} r 



drops 




Figure 12. Graphs of the relative errors E(q) (semi- logarithmic 
scale), for the initial condition computed by discretizing ( fIT]) -( IT2l) 
(dotted line) and the refined initial condition computed with the 
Gauss- Newton method (full line). The dash-dot line gives the rel- 
ative difference r(q) (in supremum norm) between last two iterates at 
the end of the Gauss-Newton iteration. 



Above a critical value q c « 0.9, the initial conditions computed as indicated 
above yield travelling waves that remain practically unchanged over very long times 
when propagating along the lattice (see figure [TBI for an example). The situation 
changes drastically below q c because the travelling waves become unstable, yielding 
an amplification at exponential rate of the initially small errors made on the initial 
condition. This situation is described in figure O for q = 7 'n /25. The travelling 
wave is destroyed by a instability at t « 360 ~ 303 T, which rapidly drives the 
system into a disordered regime. These instabilities will be described in more detail 
in the next section. 
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Figure 13. Solution of ([I]) computed over long times, for the 
initial condition determined by Gauss-Newton minimization when 
q = 37r/10 ~ 0.94. The inverse wave velocity is T ~ 1.20. Up- 
per plot : bead displacements generated by the travelling wave at 
t rj 3971 (black curve), compared with the initial condition at t = 
(dots). Lower plot : spatiotemporal evolution of the interaction forces 
V'(x n+1 (t) - x n (t)) in grey levels, for t G [3760,3980]. 



3.5. Wave instabilities. In this section we consider the dynamical equation ([I]) of 
the granular chain with N = 100 particles and periodic boundary conditions, and 
numerically study the stability of the travelling wave solutions of section 13.41 We 
recall that these travelling waves take the form (1401) . where u(.;q) is the solution 
of (I11I) - (I12I) computed numerically for a given value of q, and where we fix a = 
\\u(.;q)\\~ to renormalize the solution. We rewrite equation ([1]) in the compact 
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Figure 14. Upper plot : initial bead displacements corresponding 
to the numerical travelling wave solution with q = 77r/25 ~ 0.88. 
The inverse wave velocity is T ~ 1.19. Lower plot : spatiotemporal 
evolution of the interaction forces in grey levels. 



form X = f(X) with X(t) = (Xi(t), . . .,X N (t)Y e R 2N and X n 



[X r 



x n ), and 



denote by X^ q '(t) the solutions corresponding to travelling waves. 

Since the solution X^ g \t) is time-periodic with period r(q) = 2iia~^ , its stability 
can be analyzed using Floquet theory. Let us denote by TZ(t; to) the resolvent matrix 
of the linearized equation X = Df{X^ q '{t))X and $ g = TZ(r(q);0) the associated 
monodromy matrix. The solution X^ (t) is called spectrally stable if all eigenvalues 
of $,, lie on the unit circle, and X™({) is unstable if there exists an eigenvalue of 
modulus strictly larger than 1 (see e.g. [I0JEI]). In other words, the spectral radius 
p($ q ) determines if X^ q '(t) is unstable (p($ g ) > 1) or spectrally stable (p($ g ) = 1). 
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Since X^ q \t) is a periodic travelling wave and periodic boundary conditions are 
used, one can equivalently determine wave stability using a monodromy matrix mod- 
ulo shifts. This approach reduces the length of numerical integration significantly 
at low wavenumbers, since the linearized equations are integrated over an interval 
of length T(q) = ^T~(q). Moreover, it describes more conveniently the growth of 
perturbations in a frame moving with the waves in the linear approximation. 

To be more precise, let us denote by S the spatial shift S X = (X 2 , . . . , X N , Xxf 
associated to periodic boundary conditions. The initial condition 1™(0) corre- 
sponding to a travelling wave solution is a fixed point of the map M q defined by 

Ar q (X(0))=SX(T(q)), 

where the inverse wave velocity T{q) is given by (!60l) . Let us introduce the mon- 
odromy matrix modulo shift JF q = DJ\f q (X^ q \0)), which takes the form 

T q = SK(T(q);0). 

Recalling that q = 2mn/N (1 < m < N/2), we have r(q) = ^T(g). Thanks to 
the invariance S X {q \t + T{q)) = X®{t), one has SK{t + T; t ) = K(t\ t - T) S, 
which implies 

-pN fom 
• r q ^ q ■ 

Consequently one has p(& q ) = p{F q ) , i.e. X^(t) is unstable for p(J r g ) > 1 and 
spectrally stable if p(J r g ) = 1. 

The spectral radius of T q is plotted in figure [15] as a function of q (these results 
are obtained using the software package Scilab, for which matrix eigenvalue com- 
putations are based on the Lapack routine DGEEV). These results show that the 
travelling waves are unstable below the critical value q c ~ 0.9. The growth rate 
of the instability in the linear approximation is maximal for q = q m w 0.56, and 
decreases rapidly when q — ¥ and q — ¥ q c . Figure [16] shows the evolution of the 
eigenvalues of T q when q is varied in this parameter region. One can see that the 
number of unstable modes rapidly increases below q = q c . 

According to figure E] we have N(q c ) ~ 3, i.e. instabilities show up when the aver- 
age number N(q) of adjacent interacting beads becomes > 3. This case corresponds 
to the maximal number of adjacent interacting beads becoming > 4. This phenom- 
enon can be intuitively understood through the results of [17], where the travelling 
wave stability was numerically established for a three-ball chain (with fixed center 
of mass) using Poincare sections. However, as it follows from our numerical results, 
the interactions of the three-ball packets with additional beads generate instabilities 
as soon as q < q c . 

The weakness of the linear instability for q pa is consistent with the conver- 
gence of the travelling waves towards the classical Hertzian solitary wave (section 
13. 2p . given the fact that the Hertzian solitary wave appears structurally robust in 
dynamical simulations. 

Now let us describe how these instabilities act on the wave profiles. When q is 
slightly below q c , the early stage of the instability generates two large regions of 
average compression separated by two large regions of average extension (see figure 
IPTj) . This transitory state disappears rapidly and a disordered regime settles, as 
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shown previously in figure [HJ When q reaches q m , the large-scale transitory state is 
not visible any more and the system directly evolves from the regular travelling wave 
towards the disordered regime. The instability of the travelling waves is observed up 
to the long wave limit, where it manifests differently from the instability at q ~ q c . 
A slow dispersion of the compression pulse occurs (see figure [T8| top right plot), and 
subsequent collisions of dispersive waves with the shock generate a fast instability. 
A closer view of the transition from dispersive travelling waves to spatiotemporal 
disorder is provided by the lower plot of figure [£BJ In these different examples, 
we interpret the abrupt transition to a disordered regime as the result of the large 
number of unstable modes. 

It is interesting to note that the disordered regime that occurs after these in- 
stabilities coexists with some partial order, because transitory large-scale organized 
structures appear intermittently. This phenomenon is illustrated by figure [19] (lower 
plots). These structures seem to result from the interactions of travelling waves that 
coexist with the disordered dynamics. This is shown in the upper plot of figure [THJ 
which provides the spatiotemporal evolution of bead displacements in grey levels. 
The region corresponding to £ £ [0, 140] (nearly uniform at the scale of the figure) 
corresponds to the propagation of an almost unperturbed periodic travelling wave 
with wavenumber q = 97r/50 ~ q m , and spatiotemporal disorder occurs for t > 140 
after an instability. In this region, the presence of stripes reveals counter-propagating 
travelling waves. 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 



Figure 15. Spectral radius of the monodromy matrix modulo shifts T q . 

Above the critical value q c rj 0.9, we have observed very slow modulational in- 
stabilities by integrating ([[]), starting from small random perturbations of travelling 
waves with specific wavenumbers. However, the above Floquet analysis is not suf- 
ficiently precise to correctly account for these very small instabilities in the linear 
regime, because they may be overhelmed by numerical errors. 

4. Discussion 

Even though the Hertzian granular chain is commonly denoted as a "sonic vac- 
uum" , we have shown that long granular chains sustain periodic travelling waves on 
a full range of wavenumbers. We have proved the existence of a family of periodic 
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FIGURE 16. Eigenvalues of the monodromy matrix modulo shifts 
T q (marks), plotted in the complex plane for decreasing values of q : 
q — || (upper left plot), q — || (upper right plot), q — |^ (lower left 
plot), g = ^ (lower right plot). The unit circle is also represented. 



travelling waves with wavenumbers q close to 7r and profiles close to binary oscilla- 
tions. Using numerical continuation in q, we have been able to follow this branch 
of solutions up to the long wave limit q ~ 0, where they become close to a solitary 
wave inside small compression regions. 

The waves we have obtained display unusual properties, due to the fully nonlinear 
and unilateral character of Hertzian interactions. Each bead periodically undergoes 
a compression phase followed by a free flight, a transition associated with a limited 
smoothness of the wave profile (i.e. the corresponding solutions of (TTBl) are C 3 but 
not C 4 at the onset of free flight). Moreover, below a critical value q = q^ ps 1.8, the 
waves can be considered as a periodic train of independent compactons separated 
by beads in free flight. These numerical findings imply the existence of an isolated 
compacton in granular chains, when beads are separated by equal gaps outside the 
compression wave, the gaps depending on the velocity and width of the compacton. 

An interesting open problem is to prove analytically the existence of the periodic 
travelling waves obtained numerically, far from the limit of binary oscillations. For 
q ~ 0, this problem is equivalent to the existence of a family of compactons close to 
the Hertzian solitary wave, in the limit when the gaps between beads become small. 
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Figure 17. Bead displacements for the same travelling wave as in 
figure [HI (q = 77r/25 ~ 0.88), showing the growth of an instability at 
three different times t x ps 362, t 2 p^ 381 and t 3 ps 405. 



From a dynamical point of view, below a critical wavenumber q c ps 0.9, we have 
observed fast instabilities of the periodic travelling waves leading to a disordered 
regime. This threshold was attained when the average number of adjacent inter- 
acting beads becomes > 3, or equivalently when the maximal number of adjacent 
interacting beads becomes > 4. The travelling waves appeared far more robust 
for q > q c , i.e. when the number of adjacent interacting beads was always < 3. 
This can be intuitively understood through the results of [37], where the travelling 
wave stability was numerically established for a three-ball chain with fixed center of 
mass. However, as soon as q < q c , we have shown numerically that the interactions 
of the three-ball packets with additional beads generate instabilities. These insta- 
bilities persist up to q ~ 0, where they display a slower growth rate in the linear 
approximation. 

As illustrated by figure HH a cumulative instability effect is present due to peri- 
odic boundary conditions, since all perturbations left behind the compression pulse 
interact subsequently with it. If the existence of isolated compactons can be mathe- 
matically established, it would be interesting to analyze their spectral stability and 
determine if the above exponential instabilities persist. More generally, the spectral 
stability analysis of periodic travelling waves would deserve more investigations. 
Above the critical value q c ~ 0.9, it would be interesting to determine for which 
wavenumbers the waves are stable or slowly unstable. This question requires more 
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Figure 18. Top left plot : initial bead displacements at t — 
corresponding to the travelling wave solution with q = n/50. Top 
right plot : bead displacements at time t ~ 577 (the inverse wave 
velocity is T ~ 1.08). Lower plot : spatiotemporal evolution of the 
interaction forces in grey levels for t G [750, 1150], near the onset of 
an instability. 



refined numerical methods to resolve very slow instabilities present in this parameter 
regime. In addition, unusual perturbations of Floquet eigenvalues may arise from 
the limited smoothness of Hertzian interactions. 

Another question concerns the generation of stable periodic travelling waves in 
driven granular chains, in relation with possible experimental realizations. In prin- 
ciple, the travelling waves we have analyzed could be generated in finite systems, 
provided the motions of the first and last beads are imposed (following an exact 
travelling wave profile) and starting from an exact initial condition. However these 
conditions are extremely restrictive. In addition, dissipative effects should be also 
considered for practical applications. In this context, it would be interesting to 
study how to generate periodic travelling waves from simpler initial conditions and 
driving signals in dissipative granular systems. 
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Figure 19. Case of a highly unstable travelling wave with q = 
97r/50 ~ 0.56 (the inverse wave velocity is T ~ 1.14). Upper plot : 
spatiotemporal evolution of the bead displacements in grey levels. 
Lower plots : bead displacements at t ~ 410 (left plot) and t ~ 561 
(right plot), revealing an intermittent large-scale organized structure. 



Appendix A. Compression solitary waves 

In this appendix we recall some classical properties of solitary waves in granular 
chains which are used in section 13.21 for the analysis of the long wave regime. 

Consider travelling wave solutions of ([T]) taking the form x n (t) = y(s), where 
s = n — ct. The function y satisfies 

c 2 y"{s) = V'(y(s + 1) - y(s)) - V'(y(s) - y(s - 1)), seR, (61) 

where we recall that V'(x) = — \x\ a H (— x), a > 1 and H denotes the Heaviside 
function. Up to rescaling y, one can fix c = ±1 in 061 p without loss of generality. In 
that case, the renormalized relative displacements r(s) = y(s + |) — y(s — |) satisfy 

r"(s) = V'(r(s + l))-2V'(r(s)) + V'(r(s-l)), s6i (62) 

There exists a negative solution of (162]) satisfying lim r(s) = and r'(0) = 0, which 

s— >±oo 

is known to decay super-exponentially [211 [36], EHl HSJ HE] • This solution corresponds 
to an exact solitary wave solution of (II]) close to Nesterenko's approximate solution, 
with velocity equal to unity. Figure [20] shows the solitary wave profile computed 



32 GUILLAUME JAMES 

numerically, using the same numerical scheme as for periodic waves (except we 
change the boundary conditions and use an explicit approximation of r derived by 
Ahnert and Pikovsky [I] to initialize the Broyden method). Our results agree with 
the ones of references [H [151 EH], in particular the solution we obtain is even in s. 
To recover y from r, we note that the Poisson equation 

-y"(s) = f(s), s e R, 

admits for all odd functions / decaying exponentially at infinity a unique odd and 
bounded solution given by y(s) = L G(s, t) f(t) dt, where 

G(s,t) = -H(st){\t + s\-\t-s\). 
Moreover, one has 

r+oo 

lim y(s) = ± / tf(t)dt, (63) 

the convergence being exponential. Equation (1611) with c = 1 can be rewritten 

y"(s) = V'(r(s+ 1 -))-V'(r(s- 1 -)), 

where the right side is odd in s due to the evenness of r. Consequently, we obtain 
a solution of (IBTj) with c = 1, given by 

Vis) = [ G(s, t) [V'(r(t - -)) - V'(r(t + -))] dt, (64) 

and satisfying 

lim y(s) = k a , lim y(s) = -k a , (65) 

s— >■— oo s— >+oo 

where we have 

/■+oo 

k a = - V'(r(s)) ds (66) 

Jo 
by virtue of (163^) (this simplification is obtained using the evenness of r and ele- 
mentary changes of variables in the integral). Note that a simpler formula can be 
derived for the numerical computation of y, using the fact that 

N 1 

y{s) =y(s-N-l) + J2r(*-k- ^)- 

k=0 

Letting iV — > +oo and using the evenness of r yields 



+oo -. 

y(s) = k a + ^r(k + --s). 



k=0 

Then setting s = gives 



+oo -. 

-Y,r{k + -\ (67) 



fc=0 

and consequently 

+oo 1 1 

V(s) = ^r{k + --s)-r{k + -). (6J 



2 ' x 2' 

k=0 
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Using formula fl67|) we numerically obtain ^3/2 ~ 1.3567. 
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Figure 20. Solitary wave solution of ( 1621) computed numerically 
(upper plot) and corresponding bead displacement y(s) solution of 
( IBTj) for c = 1 (lower plot). 
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